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ABSTRACT 

Stationary solutions to the equations of non-linear diffusive shock acceleration play a fundamental 
role in the theory of cosmic-ray acceleration. Their existence usually requires that a fraction of the 
accelerated particles be allowed to escape from the system. Because the scattering mean-free-path 
is thought to be an increasing function of energy, this condition is conventionally implemented as 
an upper cut-off in energy space — particles are then permitted to escape from any part of the 
system, once their energy exceeds this limit. However, because accelerated particles are responsible 
for substantial amplification of the ambient magnetic field in a region upstream of the shock front, 
we examine an alternative approach in which particles escape over a spatial boundary. We use a 
simple iterative scheme that constructs stationary numerical solutions to the coupled kinetic and 
hydrodynamic equations. For parameters appropriate for supernova remnants, we find stationary 
solutions with efficient acceleration when the escape boundary is placed at the point where growth 
and advection of strongly driven non-resonant waves are in balance. We also present the energy 
dependence of the distribution function close to the energy where it cuts off - a diagnostic that is in 
principle accessible to observation. 

Subject headings: acceleration of particles — shock waves — methods: numerical — ISM: cosmic rays 
— ISM: supernova remnants 



1. INTRODUCTION 

It is widely thought that the diffusive acceleration of 
charged part icles at shock fronts can be very efficient (for 
a review see iMalkov fc DrTinl[200iT ). The" precise value 
of the efficiency is controlled by the microphysics of the 
injection of low-energy particles into the acceleration pro- 
cess, i.e., by processes operating on small length-scales. 
However, in the case of non-relativistic shocks, such as 
those encountered in supernova remnants, the bulk of 
the energy in nonthermal particles is carried off by those 
of the highest attainable energy, i.e., by particles that 
interact with relatively large length-scale structures. 

Conventionally, the physics of this system is captured 
by combining a hydrodynamic description of the back- 
ground plasma with the diffusion-advection equation 
obeyed by the distribution function of the accelerated 
particles. Analysis of the stationary solutions of these 
equations is the foundation on which the study of the 
overall efficiency of the process and the maximum en- 
ergy to which particles can be accelerated rests. 

The imp ortance of stationary solutions was realized 
early on by iDrurv fc Vo clk (1981) who used a reduced, 
two-fluid description to analyze acceleration by plane 
shock fronts in a one-dimensional flow. Provided accel- 
erated particles are not permitted to leave the system, 
the two-fluid description can be derived from the full 
description, including the diffusion-advection equation, 
using only two plausible assumptions about the particle 
distribution function. The more important of these is 
that the so-called effecti v e di ffusion coefficient is always 
positive. IDrurv fc Voelkl (|1981h proved that this restricts 
the possible stationary flow patterns to those containing 



a precursor in which the accelerated particles deceler- 
ate the incoming upstream flow, followed by a hydrody- 
namic shock front. This important result, which implies 
that only a finite number of stationary solutions exist for 
shocks of a given Mach number, was subsequently gen - 
eralized to the relativistic case bv lBaring fc Kirkl (|1991h . 

However, the relevance of these studies is questionable 
if the diffusion coefficient is an increasing function of par- 
ticle energy, as is the case, for example, if the trans- 
port is dominated by scattering off Alfven waves via 
the cyclotron resonance. One reason for this is that the 
timescale on which a stationary solution is approached 
becomes large at high particle energy. In this case, quasi- 
stationary solutions with a constant distribution function 
below a slowly evolving upper cut-off in energy can be ex- 
pected to es tablish them selves, and have been found nu - 
merically by [HI (|1987l) and bv lFalle fc Giddingsl (|1987l) . 
Another reason is that, as particles are accelerated to 
higher and higher energy, their mean-free path increases, 
and at some point becomes comparable to the size of any 
realistic system. Thus, even if one considers only strictly 
stationary solutions, the escape of high-energy particles 
appears to be an important property that will limit the 
maximum energy to which particles can be accelerated 
in any given system. 

One way of accounting for particle escape is to trun- 
cate the distribution function above some finite value of 
the energy. But even with this simplification, finding sta- 
tionary solutions of the combined kinetic and hydrody- 
namic equations is much more difficult than solving the 
two-fluid system. An important advance was the discov- 
ery of an approximate analytic s olution of the diffusion- 
advection equation by IMalkov) ( |l997ai ). This solution 
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has a particle distribution function that vanishes above 
an upper cut-off, p m ax, in the magnitude of the particle 
momentum. Wherever the hydrodynamic flow is com- 
pressed (usually throughout the precursor and at the 
shock front) there exists a flux of particles across this 
boundary in momentum space. These particles cease to 
contribute to the stress-energy tensor, and, therefore, ef- 
fectively escape f rom the sys t em. A similar distribution 
was adopted by lAchterbergl (|1987T ) when investigating 
numerical solutions to this equation. Under this assump- 
tion it is possible to find approximate analytic solutions, 
and relatively straightforward to const ruct numerical so- 
lutions to the full set of equations (lAchterbergl 119871 : 
iMalkov et aT]|2000t lBlas3l200a lAmato fc Blasill2005fl . 

Two problems intrinsic to this approach are immedi- 
ately obvious. Firstly, the momentum dependence of 
the distribution close to p mayi — a diagnostic that is, 
in principle, accessible to observation — is not well- 
approximated. Secondly, the exclusion of a post-cursor, 
although made plausible by the two-fluid approach, is not 
always justified. The conditions that must be fulfilled by 
the distribution function and the momentum-depe ndent 
diffus ion coefficient for this assumption to be valid (|Kirkl 
1990) can in principle be checked a posteriori, but this is 
not straightforward for discontinuous distributions. 

These problems arise because an upper cut-off in mo- 
mentum is used to describe the physics of particle escape. 
Recently, however, observational evidence has been accu- 
mulating suggesting that cosmic rays are responsible for 
substantial amplification of the ambient magnetic field 
in the precursors of shock fronts in supernova remnants 
dHwang et al.ll2002t IVink fe Laming) 120031: iBamba et all 
l2005t lUchivama et al.ll2007f ). This implies that the scat- 
tering mean-free-path is not only a function of energy, 
but also depends strongly on position with respect to 
the shock front. It, therefore, highlights a more serious 
short-coming of the approach that uses a cut-off in mo- 
mentum space: The amplification of the field is likely to 
be connected with the spatial flux of energetic particles, 
which is artificially distorted if particles are assumed to 
vanish across a momentum boundary. In this paper we 
examine stationary solutions with escape through a spa- 
tial boundary instead of through a boundary in energy 
space. 

Spatial boundaries have previously been implemented 
in Monte-Carlo simulati ons of the non-line a r acce leration 
problem. In particular, IVladimirov et al.1 (|2006| ) used a 
model equation to describe Alfvenic turbulence and cou- 
pled it to a Monte-Carlo simulation that permitted es- 
cape over a spatial boundary. In this way they were able 
to investigate the effects of an enhanced resonant inter- 
action between the turbulent waves and the accelerated 
particles, although the position of the boundary itself 
remained ar bitrary. 

Recently, iZirakashvili fc Ptuskinl (|2008fl used MHD 
simulations to describe the excited turbulence and find 
the diffusion coefficient of the highest energy accelerated 
particles. They allowed these particles to escape over a 
spatial boundary, but used the test-particle approxima- 
tion in which the flow speed is unaffected by the particles. 
In discussing shock fronts in supernova remnants, they 
suggested that the boundary leads the shock front by a 
distance given roughly by the radius of the remnant. 

Our approach differs from these, not only in the nu- 



merical method used to solve the acceleration problem, 
but also in the input physics. Resonant interactions be- 
tween energetic particles and Alfven waves were long 
thought to be responsibl e for coupl i ng these particles to 
the background plasma (iBe ll 1978; MacKeuz ie fc Voelkl 
119821 : lAchterbergl 119831 : iLucek fc Belli 120001 ) . However, 
well ahead of the shock front, non-resonant processes 
are more strongly driven and can be expected to domi- 
nate under the conditions present in supernova remnants 
(lBellll2004t|Pelletier et alj|2006t iBvkov fc Toptvginl l2005l: 
I Reville eta l. 2007; Zirakash vili et alj |2008). In the linear 
phase, these instabilities inject short-wavelength turbu- 
lence into the plasma, resulting in a relatively large dif- 
fusion coefficient that is proportional to the square o f 
the particle momentum (|Zirakashvili fc Ptuskinl [2008) . 
However, the non-linear evolution includes not only a 
cascade of energy to smaller length scales, but als o the 
appe arance of large-scale structures such as cavities (jBelll 
2005). In this regime, the mean-free-path is reduced, and 
its p dependence eliminated, as can be seen from large- 
scale numerical sim ulations of the transport properties 
(jReville et al.| [2008). These simulations have not yet ad- 
vanced to the stage where they can provide a model diffu- 
sion coefficient over a wide dynamic range of momentum. 
Consequently, we model this effect by assuming the dif- 
fusion to be of Bohm type in a background field that 
is amplified to the value at which the non-resonant in- 
stabilities are expected to saturate. Although this is a 
relatively crude approach, it enables us to solve the non- 
linear problem of finding stationary solutions. We are 
then able to check the location of the spatial boundary, 
which should be located where the growth rate of the 
non-resonant modes is approximately equal to the speed 
of advance of the shock divided by the distance from the 
shock front. Since the escaping flux is dominated by the 
highest energy particles, we do not expect that changing 
the form of the diffusion coefficient will alter significantly 
our conclusions, provided it remains an increasing func- 
tion of momentum. 

The paper is set out as follows: In Section [2] we set up 
the advection-diffusion equation and the hydrodynamic 
equations governing the system of accelerated particles 
and background plasma. The two ways of allowing for 
particle escape are discussed in Sectional In Section[4]we 
describe the iteration and finite difference methods used 
to find stationary solutions of the combined advection- 
diffusion and hydrodynamic equations using the two dif- 
ferent boundary conditions, and in Section [5] we com- 
pare and discuss the results. A discussion of the self- 
consistency of the position of the spatial boundary and 
a summary of our conclusions is presented in Section [6l 

2. BASIC EQUATIONS 

We consider a gas subshock located at x = with a 
flow profile, in the subshock rest frame, given by 

U(x) = ! T\ X < n 
v ; 1 u(x) x < 

with u(x = 0~) = u\ and 1*2 constant in the absence 
of a post-cursor. The gas velocity far upstream is de- 
noted by uo- We assume that, as a result of scatter- 
ing centres frozen into the flow, energetic particles, with 
speeds v 3> u(x), undergo diffusion with a momentum- 
dependent diffusion coefficient n(p). These particles are 
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also advected with the flow, adiabatically compressed 
and injected from the thermal background. The isotropic 
part of t heir phase sp ace density obeys the transport 
equation (Skil linelll975t ) 

— + U—-—(n—\--—p— + Q{xp) (1) 
dt dx dx \ dx J 3 da; dp ' ' 

where Qo describes the injection of particles into the ac- 
celeration process. For mono-energetic injection at the 
gas subshock 

Qo(x,p) = —^-^- S(p - p )5{x), 

where the number density of gas particles entering the 
shock front is n gi i and r\ is the fraction of entering par- 
ticles that take part in the acceleration process. In our 
notation, the particle mass and momentum are m and 
mcp, so that p, po etc. are dimensionless, and the phase- 
space distribution function / has the dimensions of an 
inverse volume. 

Integrating Eq. ([I]) first across the shock and then 
across the injection momentum, it follows that 

3ui 7?n g ,i 



h(Po) = 



Am 4ttPq 



(2) 



with Am = Mi — U2- An important restriction on this ap- 
proach is that the distribution function at the injection 
momentum po must be approximately isotropic. This 
requires that the velocity of these particles should be 
several times greater than that of the upstream plasma. 
Since we will be interested primarily in shocks in super- 
nova remnants, where Mo/c < 1/30 we require pq > 0.1. 

In this paper we are interested in steady state solutions 
to the particle transport equation, and in particular the 
role played by particle escape, when the pressure associ- 
ated with the energetic particles 



4tt 



Pct(x) = —mc / vp f(x, p) dp 



(3) 



reacts on the flow. Sufficiently far upstream, in the ab- 
sence of cosmic rays, the gas has a density po and pressure 
Pgfi. Mass and momentum conservation give 

p(x)u(x)=p u , (4) 



P CI (x) + p(x)u(x) 2 + P gfi 



u(x) 



where 7 is the adiabatic index of the gas. 

The plasma flowing towards the shock is adiabatically 
compressed and slowed in the precursor. It's velocity 
profile u{x) is a monotonic function. Non-adiabatic heat- 
ing is potentially impo rtant (e.g. ICaprioli et al.l 120081 : 
I Vladimirov et all 12008) . However, in the interests of 
simplicity, we shall neglect it in the following. For adi- 
abatic heating alone, the sub-shock compression ratio 
r s = ui /M2 and the pre-compression 

R = u /ui. (6) 

are related by (see, for example lLandau fc Lifshitzlll959f) 

7 + 1 



7 - 1 + 2Rf+ 1 M~ 



(7) 



where M is the Mach number of the flow at upstream 
infinity. For a given pre-compression ratio, we can use 



Bernoulli's equation ([5]) to determine the cosmic-ray 
pressure at the shock 



Pcr(0~ 



1 - RT 



1 



7M 2 



(1 - i? 7 ) 



(8) 



Thus, given the upstream conditions and R, or, alter- 
natively, r s , we can determine the cosmic-ray pressure at 
the shock. 

3. PARTICLE ESCAPE 
3.1. Boundary in energy 

As mentioned in Section [TJ most previous analytic and 
semi-analytic calculations adopted a precursor of infinite 
spatial extent. Escape is permitted by assuming that 
particles with energy above a certain threshold decouple 
from the plasma. In essence, this is equivalent to as- 
suming that the mean-free-path to scattering is energy 
dependent, and, above a certain threshold, becomes large 
compared to the size of the system. In the resonant scat- 
tering scenario, this implies that the relevant wave spec- 
trum cuts-off above a certain wavelength. 

Since the most energetic particles are highly relativistic 
an upper cut-off in the energy is equivalent to one in the 
magnitude of the momentum. For p > p max , particles 
escape from the system, and the approximate steady- 
state solution for p < p max is 



f(x,p) = /o(p)exp 



3k 



1 



"2 



dx u(x) 



(9) 



where q = dlnfo/dlnp. The addi t ional factor of 
(1 — M2/M1) was included in Blasi et al. (2007) to match 
the boundary condition for weakly modified shocks. The 
distribution at the shock, for monoenergetic injection, is 

, , v 3Rr s V n \ f p dp' 3Rr s U{p' 

hip) = 



Rr s U{p) - 1 Airpl 6XP 



where U(p) = u p (p)/uo with 



P 



Rr s U{p') - 1 
(10) 



Up(p) = ui - 7- / dx^-f(x,p). 
Jo J-00 dx 



(11) 



Following [Blasi (|2002f ) we take U(p ) ps R^ 1 to give 

r.-l 



(12) 



3n Rr s 

in agreement with equation In the results of this 
work, we will adopt the dimensionless in jection parame- 
ter, v, as defined in lMalkov et al.1 (|2000D . 

4?r mc 2 4 
v = ~5 2^o/o(Po), (13) 

where again, p is dimensionless. 

3.2. Boundary in space 

A free-escape boundary upstream implies that all par- 
ticles that cross a surface placed at a distance L csc up- 
stream of the shock leave the system. Essentially, the 
particle mean-free-path and, therefore, the diffusion co- 
efficient, become infinite at the escape boundary due to 
the absence of scattering waves beyond that point. Since 
we employ a diffusion approximation, this can be imple- 
mented by setting the isotropic part of the distribution to 
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zero on the boundary. Whilst inside, we assume particles 
undergo Bohm diffusion, 

k = K Q p, x > -L csc (14) 

Higher energy particles have a larger mean- free-path and, 
therefore, propagate from the shock to the boundary 
more easily than low energy ones. Since they then es- 
cape, a turn-down in the spectrum results. However, 
this is a smooth decrease, rather than an abrupt cut-off. 
It occurs close to a momentum p* defined by 

L csc = K,(p*)/u , (15) 

In the test particle limit, where the upstream flow is 
unmodified and constant, a straightforward calculation 
shows that the spectrum at the shock, fo(p), has a slope 
given by 



dln/o 3«i 



1 



dlnp Am [1 — exp(— L 6SC mo/k)] 



(16) 



At low momenta, p <C p* , the spectrum agrees with the 
standard test-particle result without escape upstream, 
but at high momenta p ~ p* the spectrum cuts off ex- 
ponentially. It should be noted that some particles with 
momenta p < p* also reach the escape boundary, and 
the spectrum will start to turn over at these values. Al- 
though the particle distribution function / vanishes at 
the free-escape boundary, the particle flux remains finite, 
being proportional to p 2 ndf / dx. In the test-particle ap- 
proximation it can be written, making the substitution 
s =p/p*, 



S 2 K 



df _ u Q fo(p )s 2 
dx ~ e 1 /* - 1 



■ exp 



-3r s 



1 



din; 



1 



-l/s' 



17) 



where s = po/p*. As a function of momentum, the 
escaping flux is sharply peaked atsw (r s — l)/(f s + 2). 

4. METHOD OF SOLUTION 

In the case of a boundary in energy, we use an itera- 
tion s cheme similar to that employed bv lAmato fc Blasl 
(2005). Given a shock Mach number M, shock speed uq 
and momentum range po, p ma x, we look for converged 
solutions for each value of R. The initial flow profile is 
linear, and the spectrum is a single power law. This al- 
lows us to calculate U(p) from which we determine the 
injection parameter rj in Eq (jlOj) . by identifying the cos- 
mic ray pressure at the shock front in Eq ([5]) with the 
integral over p of vp 3 fo, as indicated by Eq ([3]). Using 
equations (J9j) and ((3]) we then update the flow using the 
fluid equations. This gives new values for u(x) and q{p), 
and, hence a new value for r/. This is repeated until r/ 
has converged. 

Numerical solutions for steady state modified shocks 
and particle spectra, with a free escape boundary up- 
stream, are also found using an iterative procedure. For 
a given Mach number (M), shock speed (uq) and pre- 
compression ratio (i?), we initialize the spectrum using 
the test particle solution. The flow profile is then ad- 
justed using the flux conservation equations (|4I5|) . With 
the modified flow, the distribution function is updated by 
solving the time dependent transport equation for par- 
ticles, equation |T]), with the upstream boundary condi- 
tion f(—L esc ,p) — 0. For this we use a Crank-Nicholson 
scheme centred in space and upwind in momentum. Us- 
ing Eq. (|8l), the distribution function can be normalized 




Fig. 1. — The injection parameter u as a function of the pre- 
compression R — see Eq. — for increasing Mach number at 
a fixed upstream velocity no = 5, 000 km s~ 1 . In generating this 
plot, we set po = O.lmc and p* = p max = 10 3 mc. The solid lines 
are found using a spatial boundary, the dashed lines use a boundary 
in energy. 



to match the pressure at the shock. When, at each iter- 
ative step, the transport equation is solved, the resulting 
cosmic ray pressure is used to update the flow which, in 
turn, is then used to update the particle distribution. 

The solution is found when the fluid quantities no 
longer change and the injection parameter r\ has con- 
verged. The code steps through values of the pre- 
compression ratios, using the previous converged profile 
and spectrum as the initial conditions for the next iter- 
ation. In this manner we obtain a full set of solutions 
that depend on rj, M, R and L esc - 



5.1. 



5. RESULTS 

Comparison of approaches 

In order to make a meaningful comparison between the 
two methods we take 



P' =Pn 



(18) 



so that the boundary for upstream spatial escape, char- 
acterized by L csc , is placed one diffusion length away 
from the shock for particles with momentum p ma x- Con- 
sequently, the gradual momentum cut-off produced by 
the spatial boundary technique is close to the position of 
the abrupt momentum boundary. 

Figure[T]plots the numerically determined injection pa- 
rameter v as a function of the pre-compression ratio for 
the two different methods. These plot s follow a simila r 
form to those studied analytically by iMalkovl (|1997bf) . 
The curve for each shock Mach number can, in general, 
be separated into three different regimes. The inefficient 
weakly modified regime where the spectrum resembles 
the test particle solution, the efficient regime where the 
shock is strongly modified, with a weak subshock, and 
an intermediate regime between the two. This inter- 
mediate regime lies typically in a narrow range of the 
injection parameter v. As the shock Mach number in- 
creases, mul tiple solutions fo r a single valu e of v occur, 
as found bv IMalkovl (|1997bD ; lAmato et all (|2008l ). The 
two approaches give similar results, as expected, since 
the normalization of the spectrum is fixed by Eq ([8]). 
The differences arise due to the shape of the spectrum 
close to the cut-off. In Fig [2] we show the v — R diagram 
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R 

Fig. 2. — The injection parameter v as a function of the pre- 
compression ratio R for a M = 50 shock using different injection 
and escape momenta.. Solid: po = O.lmc, p* = 10 3 mc, dashed: 
Po = O.lmc, p* = 10 4 mc. 



for a Mach 50 shock for different values of maximum and 
minimum momentum using the free-escape boundary. 

Figures [3] and [?] compare the flow and pressure pro- 
files, as well as the spectrum and spectral index for Mach 
numbers, M = 100 and M = 500. Both examples are 
with a pre-compression ratio of R = 20. As can be 
seen from Fig. [TJ the two cases fall in different regimes. 
For the M — 100 case, the shock is strongly modified 
with r s w 2.1, while the M = 500 case is in the in- 
termediate regime with a subshock compression ratio of 
r s rs 3.86. This can be seen from the low energy shape of 
the spectrum: the reduced sub-shock compression leads 
to a much softer spectrum at the M = 100 shock. 

The flow profiles of the two models differ only slightly. 
This follows naturally from the boundary conditions: 
with a free escape boundary the precursor length is fixed 
by i C sc- Using a boundary in momentum leads to a 
slightly extended precursor, due to the larger number 
of particles with momenta p* interacting on a slightly 
larger scale. For both Mach numbers, there is a smooth 
turnover in the spectrum at momenta close to p* . Fur- 
thermore, the momentum derivative of the distribution 
remain s negative in this range. In this case, it can be 
shown |l\iik|[TT) !)(). Eq (12)) that the dynamics of parti- 
cles close to p* do not permit a post-cursor. However, 
the possibility remains that the dynamics of th e injection 
mechanism may still do so (jZank et al.lll993f ). Particles 
with momenta slightly less than p* are also able to dif- 
fuse a distance L osc from the shock and escape upstream. 
This is what leads to the observed turn-over. The re- 
duced value for p* is due to the non-linear effects of the 
system, and a crude estimate of the reduction is given by 
the formula 

Pc* ff = ^T [ u&x. (19) 
uqL J l 



5.2. 



Location of the spatial boundary 

In our approach, we assume turbulence is generated 
at the spatial boundary by a non-resonant instability 
driven by particles escaping into the undisturbed up- 
stream medium. The condition that these waves are 
strongly driven is 



c 



]crP cS mc 
ep ul 



(21) 



and Ma is the Alfven Mach number of the shock, in the 
medium upstream of the boundary. 

The cosmic-ray current is evaluated from integration 
of the diffusive flux over momentum space 



3ct{x) = -Aire 



df 2a 
K^r-p dp- 
OX 



(22) 



The integrand in this function peaks close to p — p* and 
can be easily extracted from the numerical solution to 
the full non-linear problem. Because the magnetic field 
ahead of the boundary is not specified in the solutions, we 
plot (M 2 rather than £M A as a function of the precursor 
compression R in Fig. [5j In the medium surrounding a 
supernova remnant, we expect Ma ~ M. Therefore, ac- 
cording to this figure, the non-resonant waves are indeed 
strongly driven, as defined by Eq. (|20[) . In agreement 
with the linear theory, we find that j cr (— L) is an increas- 
ing function of the shock's Mach number and a decreas- 
ing function of p* . For fixed maximum momentum, this 
result is independent of the the diffusion coefficient in 
the precursor, and, therefore, of the strength of the am- 
plified magnetic field — a weaker field leads to a larger 
precursor but does not change the escaping particle flux. 

The assumption that particles escape freely upstream 
of the boundary implicitly assumes that the instability 
responsible for the generation of the turbulence operates 
on a length scale that is short compared to the precursor 
length. This requires that the inverse of the maximum 
growth rate 7 ma x of the non-resonant instability should 
be less than the advection time through the precursor. 
The maximum growth rate is related to the driving pa- 
rameter by 

_ CMauq ,„„,. 
7max "2r g fe) {I6) 

where the gyro-radius of escaping particles is evaluated in 
the ambient magnetic field Bq upstream of the boundary. 
It is convenient to define the advection length 

u 



Ladv — ' 



7n 



(24) 



However, the length scale in our numerical simulation is 
set by the value of the diffusion coefficient in the am- 
plified magnetic field B s . Relating this to the diffusion 
coefficient by assuming Bohm diffusion gives the follow- 
ing lengthscale for our hydrodynamic precursor: 



cr g (p* 
3wo 



Bo 



(25) 



The ratio of these two length scales is an important pa- 
rameter 

Lady _6VA ( Bs 

' "c c VA), 



Pcff 

p* 



(26) 



(20) 



and, in order to evaluate it, we must e stimat e the 
strength of the am plified field. Ed] (|2004h and 
iPelletier et alj (|2006f ) give the following approximation 
for the saturated magnetic field energy density: 

S»fPc«(0-). (27) 
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Fig. 3. — Normalized cosmic-ray pressure and flow profile for a M = 100 shock, as a function of distance upstream of the sub shock. The 
pre-cursor compression ratio is R = 20. The solid lines correspond to the spatial boundary approach and dashed lines to the momentum 
boundary. 




Fig. 4. — Normalized cosmic-ray pressure and flow profile for a M = 500 shock, as a function of distance upstream of the sub shock. The 
pre-cursor compression ratio is R = 20. The solid lines correspond to the spatial boundary approach and dashed lines to the momentum 
boundary. 
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Numerical investigations of the non-linear behavior of 
magnetic turbulence in the presen ce of streaming cos- 
mic rays have been performed (e.g. iNiemiec et al.ll200cl 
iRiquelme fc Spitkovsky||2008[ ). but the results concern- 
ing the saturation level of this instability are inconclu- 
sive. Adopting Eq. (|27|) , and substituting into Eq. (j26|) we 
find 

L pc ~ C V c ) { Rpoul ) \p* J ' 1 J 

We plot iadv/^pc determined numerically according 
to Eq. (|28p as a function of the injection parameter v in 
Fig. [6l for M = 100. For weakly modified shocks with 
r s 4, the non-resonant instability grows slowly ahead 
of the escape boundary: L ac \ v 3> L pc . For intermediate 
strength modified shocks 4 < R < 10 (see Fig.[T]to relate 
v to R), the growth time is comparable to the advection 
time. For highly modified shocks, R > 10, the instability 
grows very rapidly in a small region just ahead of the 
escape boundary. 

Formally, the calculations we present are valid only in 
this latter case, where the region in which the magnetic 
field increases from its ambient strength to its strength 
in the precursor is small compared to the length of the 
precursor itself. However, the point at which amplifica- 
tion sets in is arbitrary in this case. The flux of particles 
escaping from the shock front remains constant in planar 
geometry, so that we should expect the pre-cursor length 
to increase as these particles penetrate further and fur- 
ther upstream. This leads to a higher p* and reduces 
the flux of escaping particles. In a fully self-consistent 
picture, this process should regulate itself such that the 
advection length becomes comparable to the precursor 
length. 

On the other hand, for weakly modified shocks, where 
Ludv S> L pc , the escaping particles must penetrate a 
large distance in front of the shock before the instability 
they drive has time to grow appr eciably. In essence, this 
is the situation investigated by Zirakas hvili fc Ptuskinl 
( 2008), who assumed the precursor to the shock — which 
is the region where the hydrodynamics are influenced by 
the accelerated particles — was very short compared to 
the other length scales in the problem. According to our 
results, this situation is consistent with stationarity of 
the accelerated particle spectrum only for weakly mod- 
ified shocks that are relatively inefficient. In this case, 
however, the magnetic field is amplified very gradually, 
and the approximation that the amplification region is 
short compared to the pre-cursor is inconsistent. This is 
because a level of turbulence close to that assumed in the 
precursor is already available and able to interact with 
particles far ahead of the escape boundary. This again 
suggests that, in a self-consistent picture, the precursor 
length will be comparable to the advection length of the 
instability. 

According to Fig. O the value of v at which the ratio of 
^adv to L pc is unity is a decreasing monotonic function of 
p* , suggesting the physical picture evolves in the follow- 
ing manner: For a given injection parameter v and mo- 
mentum p* , determined from the length of the precursor, 
we can calculate the ratio of the advection length to the 
precursor length. If this ratio is large the escape bound- 
ary is too far from the shock and p* will reduce itself 
until the ratio reaches unity. If, however, the ratio is too 
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Fig. 5. — The driving parameter £ times the square of the sonic 
Mach number M as a function of subshock compression for different 
parameters. Over the entire parameter range f M 2 ^> 1. Therefore, 
provided the plasma /3 = M 2 /M^ is not too large, the non-resonant 
mode is, according to Eq. 1201 . the fastest growing mode. The lines 
correspond to M = 500, p* = 10 4 (dashed), M = 100, p* = 10 5 
(solid), M = 100, p* = 10 4 (dash-dot), M = 50, p* = 10 4 (dotted), 
with an injection momentum of po = 0.1. 
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Fig. 6. — The ratio of advection length to precursor length, 
^adv/^pc, for different values of p* for a shock Mach number of 
100. From left to right the lines correspond to p* = 10 5 , p* = 10 4 , 
p* = 10 3 , p* = 10 2 . All models were calculated using an injection 
momentum of po = 0.1. The horizontal line corresponds to L a dv = 

small, p* will increase, as described above. In general, in 
the range of precursor compression rati os that interest us 
here, v is a decreasing function of p* ( Malkov fc Drurvl 
[200lh . which can be clearly seen in Fig. [21 It is natural to 
expect that the system will organize itself such that the 
advection length approaches the precursor length. For 
the parameters adopted in Fig. [6l this occurs for shocks 
in the intermediate range of modification. 

This scenario suggests there exists a relationship be- 
tween the injection and the maximum momentum in the 
system. A s imilar connect i on ha s previously been in- 
vestigated by iMalkov et al.l (|2000[ ) in the context of self- 
organized criticality in cosmic-ray modified shocks, al- 
though the value of the maximum momentum is con- 
trolled, in our case, by a different mechanism. However, 
as in their case, the actual solution depends on the micro- 
physics at the subshock, which determines the injection. 

6. CONCLUSIONS 

The assumption underlying most previous investiga- 
tions of non-linear diffusive shock acceleration is the non- 
existence of waves able to scatter particles with energy 
above an upper cut-off. In this paper, we examine an al- 
ternative picture, in which the current generated by the 
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streaming cosmic rays falls below some critical value at a 
distance L csc upstream of the subshock, and that beyond 
this the turbulence is insufficient to scatter the particles. 

In terms of the non-linear response of the system to 
changes in the injection parameter, we find the two pic- 
tures are quite similar (see Fig.[lJ. One difference is that 
our spatial boundary method can be used to model the 
shape of the distribution close to the cut-off (see Fig. [3] 
and Fig. 2]). This is an advantage because it potentially 
enables one to model the radiative signatures of the ac- 
celeration process. 

The main difference, however, is that we are able 
to address the physics that determines the location of 
the boun dary. Previous wo r k tha t implemented such a 
boundary IVladimirov et all (|2006f ) did not constrain its 
location. In time-d ependent models o f acce leration in 
supernova remnants iBerezhko fc Vol k (1997), an effec- 
tive spatial boundary is imposed by the spherical geome- 
try. But they assume a value for the amplified magnetic 
field in the entire computational box, without consider- 
ing whether or not this is consistent with the location of 
the boundary. 

We argue that the location of the boundary is deter- 
mined by the growth rate of the instability responsi- 
ble for field amplification. It has recently been shown 
that, in the case of efficient sh ock a c celera tion, a short- 
wavelength non- resonant mode [Bell ( 200 H) play s a cru - 
cial role. In a s e ries o f papers, " iPelletier et all (|2006fl ; 
iMarcowith et al.l (|2006f ) argue that this non-resonant 
mode dominates far from the shock, while the reso- 
nant streaming instability takes over closer to the shock, 
driving the diffusion towards Bohm-type in the precur- 



sor. According to these arguments, the position of the 
free-escape boundary should be fixed by the properties 
of the non-resonant mode, once the shock is modified 
by the cosmic-ray pressure. We show explicitly that 
the non-resonant modes are strongly driven at the es- 
cape boundary, and compare the local growth rate with 
the rate at which the waves are advected towards the 
shock front (see Fig. [6J. In this picture, particles stream 
freely ahead (upstream) of the boundary, whereas be- 
hind (downstream of) it, the turbulence is assumed to 
be fully developed and particles undergo Bohm diffusion 
i n an amplified magnetic field. 

Zirakash vili &: Ptuskinl (|2008| ) also suggest that the ge- 
ometry of a supernova shock front sets the length scale 
for the precursor, and thus determines the maximum par- 
ticle energy. However, we find that when the non-linear 
dynamics of the acceleration process are included, the 
length scale is set by the strength of the amplified mag- 
netic field and the efficiency of the injection process. Ul- 
timately, it is the microphysics of injection — represented 
in Figs. [T] and [5] by the parameter v — that determines 
not only the efficiency of the acceleration process, but 
also the maximum attainable particle energy. 
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